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Abstract 

As the continuum limit is approached, lattice QCD simulations tend to get trapped in the 
topological charge sectors of field space and may consequently give biased results in practice. 
We propose to bypass this problem by imposing open (Neumann) boundary conditions on 
the gauge field in the time direction. The topological charge can then fiow in and out of 
the lattice, while many properties of the theory (the hadron spectrum, for example) are 
not affected. Extensive simulations of the SU(3) gauge theory, using the HMC and the 
closely related SMD algorithm, confirm the absence of topology barriers if these boundary 
conditions are chosen. Moreover, the calculated autocorrelation times are found to scale 
approximately like the square of the inverse lattice spacing, thus supporting the conjecture 
that the HMC algorithm is in the universality class of the Langevin equation. 



1. Introduction 

One of the central goals in numerical lattice QCD is the computation of the proper- 
ties of the hght mesons and baryons with controlled errors. While the most impor- 
tant systematic errors in these calculations (finite volume and lattice spacing effects) 
are theoretically well understood, the relevant time scales in QCD simulations re- 
main unpredictable. In practice, the correctness of the simulations within the quoted 
statistical errors can therefore only be established through empirical tests and thus 
only to a limited level of confidence. 

In order to preserve the translation symmetry, the lattice theory is usually set up 
with periodic boundary conditions in all space-time directions. A side-effect of this 
choice of boundary conditions is the emergence of disconnected topological sectors in 
the continuum limit. On the lattice the sectors are not strictly separated from each 
other, but the relative weight in the functional integral of the gauge fields "between 
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the sectors" decreases with a high power of the lattice spacing [1] . As a consequence, 
transitions from one sector to another tend to be suppressed in the simulations and 
may eventually become so rare that a proper sampling of the sectors would require 
far longer runs than are practically feasible [2-4]. 

In this paper we address both issues, the very long autocorrelation times caused 
by the emergence of the topological sectors and the lack of theoretical control over 
the simulations. The first of them we propose to avoid by imposing open boundary 
conditions on the gauge field in the time direction (see sect. 2). With this choice, the 
field space in the continuum theory becomes connected and the topological charge 
can smoothly flow in and out of space-time through its boundaries. All statistically 
relevant parts of field space are therefore expected be accessible to the simulation 
algorithms without having to cross higher and higher topology barriers as the lattice 
spacing is reduced. 

When properly renormalized, some algorithms may even converge to a well-defined 
stochastic process in the continuum limit (sect. 3). In asymptotically free theories, 
such algorithms have a predictable scaling behaviour as a function of the lattice 
spacing and are thus theoretically controlled to some extent. The HMC algorithm [5] 
recently turned out to be non-renormalizable in perturbation theory and is therefore 
not of this kind [6] , but the algorithm may conceivably fall in the universality class of 
the Langevin equation (whose renormalizablity was established long ago [7,8]). The 
empirical tests reported in sections 4 and 5 partly serve to verify that the topology 
barriers are indeed absent if open boundary conditions are imposed and partly to 
find out whether the HMC algorithm scales like an algorithm that integrates the 
Langevin equation. 



2. QCD with open boundary conditions 

Open boundary conditions are easily imposed in QCD and do not give rise to impor- 
tant theoretical complications. While the discussion in this section is more generally 
valid, the gauge group is taken to be SU(3) from the beginning and we assume there 
is a multiplet of quarks in the fundamental representation of SU(3). Our notational 
conventions are summarized in appendix A. 

2.1 Boundary conditions in the continuum theory 

In the continuum limit, the gauge and quark fields live on a four-dimensional space- 
time Ai with Euclidean metric, time extent T and spatial size Lx Lx L. Time thus 
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runs from to T, while space is taken to be a three-dimensional torus, i.e. all fields 
are required to satisfy periodic boundary conditions in the space directions. At time 
and T, the boundary conditions imposed on the gauge potential Afj_{x) are 

Fok{x)\^^^^ = Fokix)\^^^T. = for ah k = 1,2,3, (2.1) 

where F^i,{x) denotes the gauge- field tensor. Note that these conditions preserve the 
gauge symmetry and therefore do not constrain the gauge degrees of freedom of the 
field. In perturbation theory, the boundary conditions on the latter instead derive 
from the gauge-fixing procedure. If the usual Lorentz-covariant gauge is chosen, for 
example, the time and space components of the gauge potential are found to satisfy 
Dirichlet and Neumann boundary conditions, respectively!. 

In the case of the quark and antiquark fields tp{x) and tp{x), we require that 

P+i^{x)l,=, = P-i^{x)l,=T = 0, P± = W^ 7o), (2.2) 

V^W^-Lo=o=V^(^)^+Lo=T = 0- (2-3) 

These boundary conditions are familiar from the discussion of the QCD Schrodinger 
functional [9,10]. Many of the theoretical results obtained in that context can actu- 
ally be reused here. In particular, as explained in ref. [11], one is practically forced 
to choose the boundary conditions (2. 2), (2. 3) if parity and the time reflection sym- 
metry are to be preserved. The action of the theory (without gauge fixing terms) is 
then given as usual by 

S = -^f A^xii{F^,{x)F^,{x)}+ f d^xi^{x){ji,D^ + Mo)i^{x), (2.4) 

where go and Mq are the bare coupling and quark mass matrix. 
2.2 Topology of the classical field space 

Since M is contractible to a three-dimensional torus, all SU(3) principal bundles over 
M are trivializable. Smooth classical gauge potentials may therefore be assumed to 

f As is the case with periodic boundary conditions, perturbation theory in finite volume is compli- 
cated by the presence of non-trivial gauge-field configurations with vanishing action (the constant 
Abelian fields). The remarks made here and below on perturbation theory refer to the situation at 
L = oo and finite T, where the minimum of the action is unique up to gauge transformations. 
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be globally defined difFerentiable fields. In view of the non-linearity of the boundary 
conditions (2.1), the classical field space is however not a linear space. 

We now show that the field space is connected. Evidently, any given gauge poten- 
tial An{x) satisfying Ao(x) = doAk{x) = at xq = and xq = T can be smoothly 
contracted to zero, without violating the boundary conditions, by multiplication 
with a scale factor. These fields are therefore continuously connected to the classi- 
cal vacuum configuration. On the other hand, if one starts from an arbitrary field 
Af^{x) in the classical field space, a smooth curve of gauge transformations As{x), 
< s < 1, may be defined through the differential equation 

{do + sAoix)) A«(x)-^ = 0, A.(x)U„=o = 1- (2-5) 

When applied to the potential A^{x), the transformation generates a curve in field 
space (parametrized by s) along which the field is continuously deformed to another 
field at s = 1 with vanishing time component. The transformed field can then be 
contracted to zero, as explained above, which proves that the field space is connected. 

The absence of disconnected topological sectors goes along with the fact that the 
topological charge 

Q = J^d^x q{x), q{x) = e^^p^tr{Fi^^{x)Fp^{x)}, (2.6) 

is not quantized. When an instanton on Ai is contracted to the vacuum configura- 
tion, for example, the charge fiows away through the boundaries and Q smoothly 
varies from 1 to 0. It may be worth noting here that the massless Dirac operator has 
no zero modes in the space of complex quark fields satisfying the boundary conditions 
(2.2). Moreover, the eigenvalues A of the Hermitian Dirac operator 75 (7^-D^ -f- m) 
are all in the range |A| > m (see ref. [11], sect. 2.2, for a proof of these statements). 
As long as the quark masses are non-negative, the quark determinant has therefore 
a definite sign and never passes through zero even if some masses vanish. 

2. 3 Renormalization and stability of the boundary conditions 

The renormalization of quantum field theories on space-time manifolds with bound- 
aries in general requires the usual (bulk) counterterms to be added to the action as 
well as further counterterms that are localized at the boundaries [12]. In the present 
case, the symmetries of the theory, power counting and the fact that ijnl^ vanishes at 
the boundaries however exclude such boundary counterterms. The renormalization 
of the theory thus proceeds as in infinite volume by renormalizing the coupling, the 
quark masses and the fields in the correlation functions considered. 
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Boundary conditions are subject to renormalization too and sometimes require a 
fine-tuning of boundary countcrterms. Neumann boundary conditions in scalar field 
theories, for example, arc known to be unstable under quantum fluctuations [12]. 
The situation in QCD is safe from this point of view, because there are no relevant 
or marginal boundary counterterms with the required symmetries. In particular, the 
boundary conditions (2.1)-(2.3) are stable under quantum fluctuations (see ref. [11], 
sect. 3, for a broader discussion of the subject). 

2.4 Lattice formulation 

The lattice theory is set up on a hypercubic lattice with spacing a, time-like extent 
T + a and spatial size Lx Lx L, where T and L are integer multiples of a. Periodic 

boundary conditions are imposed on all fields in the space directions, while time runs 
from to T inclusive, the terminal timc-sliccs being the boundaries of the lattice. 

As usual the gauge and quark fields reside on the links and points of the lattice. 
In particular, the link variables U{x,ij) G SU(3) live on all links with both 

endpoints in the specified range of time. The Wilson gauge action is then given by 
the sum 

5G = ^^^p)tr{l-C/(p)} (2.7) 
% p 

over all oriented plaquettes p on the lattice, U (p) being the ordered product of the 
link variables around p. Only those plaquettes are included in the sum whose corners 
are in the time interval [0, T]. The weight w{p) is equal to 1 except for the spatial 
plaquettes at time and T, which have weight i. 

In the case of the fermion fields, a possible choice of the action is 

T-a 

= «^ V'(^) (^w + Mo) V(x), (2.8) 

xo=a X 

where 

= |7m (V* + V^) - |aV*V^ (2.9) 

denotes the Wilson-Dirac operator and the fields are assumed to satisfy the bound- 
ary conditions (2. 2), (2. 3). Since the action (2.8) depends on the quark fields at times 
< xq < T only, one may then just as well set all components of the fields at time 
and T to zero. The dynamical components of the quark fields are thus those residing 
in the interior of the lattice. 
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The functional integral and the basic correlation functions are now defined in the 
standard manner. Evidently, only the dynamical components of the fermion fields 
are integrated over in the functional integral. Note that the quark determinant is a 
product of real factors, one for each quark flavour, since the Wilson-Dirac operator is 
75-hermitian with the chosen boundary conditions. The established QCD simulation 
algorithms can therefore be applied straightforwardly. 

It may not be completely obvious at this point that the fields satisfy the boundary 
conditions (2.1)-(2.3) in the continuum limit. As already mentioned in subsect. 2.3, 
these boundary conditions are stable under quantum fluctuations, i.e. it suffices to 
check that they emerge at tree-level of perturbation theory when the lattice spacing 
is sent to zero. The explicit expression for the quark propagator obtained in ref. [13] 
and a similar computation of the gluon propagator in the standard covariant gauge 
actually show this to be so. 

2.5 Quantum-mechanical representation 

The formulation of the lattice theory described above admits a quantum mechanical 

description in terms of a Hilbert space Ti of physical states and a bounded, positive- 
definite transfer matrix T [14] . In particular, the partition function of the theory is 
given by the expectation value 

Z={n,T^/''n) (2.10) 

of a power of the transfer matrix in a state O G H that encodes the chosen boundary 
conditions at time and T. 

A relatively direct way of introducing the transfer matrix formalism starts from a 
representation of the physical states through wave functions depending on a gauge 
field V{x, k) and the components 

X-(x) = P-x{S), x+{x) = X{S)P+ (2-11) 

of a quark field on the spatial lattice (see ref. [10], for example). The boundary state 
Q is represented by the wave function 

n{V,X-,X+) = {det(l + aMo- iaMVfe)}"' (2.12) 

in this language, the covariant derivatives being evaluated in presence of the gauge 
field V. Note that this expression is manifestly invariant under the gauge symmetry, 
the lattice symmetries and the (vector-like) flavour transformations. In other words, 
Q has the quantum numbers of the vacuum state. 
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Correlation functions of gauge-invariant fields have a quantum mechanical inter- 
pretation as well. The two-point function of a local scalar field ^(x), for example, is 
given by 

{ct){x)ct){y)) = -(Q,T(^-^°)/"<^(x)T(^''-2"')/«0(y)T2'°/"J7), (2.13) 
xo>yo Z 

where denotes the operator field associated to 4){x) [14]. Since the transfer 
matrix and the space of physical states are independent of the boundary conditions 
at time and T, the hadron masses and many other physical quantities can, in 
principle, be extracted from such correlation functions in much the same way as on 
lattices with periodic boundary conditions in time. 

2.6 On- shell 0(a) improvement 

The 0(a) improvement of the lattice theory follows the lines of refs. [15,16]. There is 
actually very little difference with respect to the case of the Schrodinger functional 
discussed in the second of these papers. In particular, all bulk 0(a) counterterms and 
their coefficients are exactly the same as those required for the on-shell improvement 
of the theory on the infinite lattice. 

We wish to emphasize at this point that a further improvement is not needed if 
one is exclusively interested in the correlation functions of fields localized far away 
from the boundaries of the lattice, where the effects of the latter are exponentially 
suppressed. The improvement of correlation functions involving fields close to or at 
the boundaries however requires the addition of the 0(a) boundary counterterms 

55G,b = i(cG - l)i;tr{l - U{p,)}, (2.14) 



SSp,^ = (Cp - l)a=^E{^W^WLo=a + ^(^)^(^)Lo=T-a}- (^.IS) 

X 

In these equations, Ps runs over all space-like oriented plaquettes at the boundaries 
of the lattice and the coefficients 

cx = 1 + cl%2 + ci'^^o' + • • • (2.16) 

must be adjusted so as to cancel the boundary effects of order a (since the boundary 
conditions are not the same, there is no reason to expect these coefficients to coincide 
with those needed to improve the Schrodinger functional). 
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2. 7 Other lattice formulations of the theory 



Lattice QCD with open boundary conditions can be set up in many different ways. 
Universality actually suggests that the details of the lattice theory become irrelevant 

in the continuum limit if the gauge, space-time and flavour symmetries are respected. 
Lattice formulations that preserve chiral symmetry away from the boundaries exist 
as well, but some care is required in this case in order to guarantee the locality of 
the lattice Dirac operator near the boundaries [17,11]. 



3. Dynamical properties of QCD simulations 

The interpretation of simulation data requires good control over the simulation dy- 
namics. In this section, the relevant notions are briefly discussed and some specific 
issues are addressed, which arise when studying the scaling behaviour of QCD sim- 
ulations. 

3.1 Autocorrelations 

QCD simulation algorithms produce random sequences of gauge-field configurations 
recursively, where the next configuration is obtained from the current one according 
to some transition probability. The simulation algorithms considered in this paper 
are the HMC algorithm [5] and the closely related SMD (stochastic molecular dy- 
namics, or generalized HMC) algorithm [18]. In both cases the simulation time is 
proportional to the molecular-dynamics time in lattice units and will, for simplicity, 
be identified with the latter in the following. 

Let Oi be a set of real-valued unbiased observables labeled by an index i. Their 
values Oi{t) measured at simulation time t are statistically correlated to some extent, 
i.e. the connected parts of the n-point autocorrelation functions 



in general do not vanish. In this equation, the bracket ((.••)) stands for the average 
over infinitely many statistically independent parallel simulations, which is the same 
as the average over time translations if the simulation is ergodic. 

The connected parts of the autocorrelation functions tend to fall off exponentially 
at large separations in simulation time. In particular, the two-point autocorrelation 
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functions 



T,,{t) = {{os)OM)) - {{o^{t))){{o,m) 



(3.2) 



can be shown to have a spectral decomposition of the formf 



oo 



rii(i) = Re{ci„c,„Alfl }, |A„| = e-i/"' 



71 



(3.3) 



n=0 



where tq > ti > . . . are the so-called exponential autocorrelation times of the algo- 
rithm. While these are independent of the observables considered, the coefficients 
Cin measure how strongly the observables Oj couple to the eigenmode number n of 
the transition probability. Note that neither the spectral values A„ nor the coeffi- 
cients Cin are guaranteed to be real, except in the case of the HMC algorithm and 
the Langevin limit of the SMD algorithm. 

In practice, the integrated autocorrelation times 



of the observables of interest play an important role, where At is the separation in 
simulation time of the observable measurements. The sum in eq. (3.4) amounts to 
a numerical integration of the normalized autocorrelation function pi{t) using the 
trapezoidal rule. In particular, in the Langevin limit of the SMD algorithm or if the 
HMC algorithm is used, the formula 



and thus the bound T\at{Oi) < tq hold up to integration errors. 
3.2 Topology-changing transitions 

On lattices with periodic boundary conditions, the probability per unit simulation 
time for a HMC or an SMD trajectory to pass from one topological sector to another 
is a rapidly decreasing function of the lattice spacing [2-4] . Such topology-tunneling 

f The HMC and the SMD algorithm both evolve the gauge field U {x, fj,) together with its canonical 
momentum tt{x, /i) (cf. subsect. 4.1). Equation (3.3) is partly a consequence of detailed balance in 
phase space and only holds for observables that do not depend on the momentum. 




oo 



(3.4) 



Tint(Oi) 




n 



(3.5) 
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transitions are non-perturbative lattice artifacts that may informally be described 
as "an instanton falling through the lattice" . The integrated autocorrelation time 
of the topological charge consequently tends to become very large, sometimes to the 
extent that the correctness of the simulation is compromised. 

With open boundary conditions, the situation is different, because the topological 
charge can change smoothly along a molecular-dynamics trajectory by flowing in 
and out of the lattice via its boundaries. A catastrophic slowdown of the algorithms 
as in the case of periodic boundary conditions is therefore not expected. 

3.3 Renormalizable algorithms 

The n-point autocorrelation functions of gauge-invariant local fields formally look 
like the correlation functions in a field theory in five dimensions, where the simulation 

time is the fifth space-time coordinate. When the lattice spacing is taken to zero, 
the autocorrelation functions may then conceivably have a continuum limit, provided 
the fields and the parameters (of both the theory and the simulation algorithm) are 
properly renormalized. 

Algorithms that integrate the Langevin equation are known to be renormalizable 
in this sense to all orders of perturbation theory [7,8]. An example of an algorithm of 
this kind is provided by the SMD-y algorithm (cf. subsect. 4.2). The simulation time 
has physical dimension [length]^ in this case and must be renormalized according to 

t = Ztt^lo? (3.6) 

where t is the simulation time in lattice units, Zf^go) a renormalization constant and 
in, the renormalized simulation time in some physical units. Further renormalization 
is not required apart from the usual field and parameter renormalization. 

Beyond perturbation theory, the renormalizability of an algorithm (and thus the 
associated scaling laws) can break down as a result of non-perturbative lattice arti- 
facts. On lattices with periodic boundary conditions, topology-changing transitions 
have this effect in the case of the SMD-y algorithm. However, if open boundary con- 
ditions are chosen, there is currently no reason to expect that the renormalizability 
of the algorithm does not extend to the non-perturbative level. 

3.4 Scaling behaviour of the HMC algorithm 

Free-field studies of the HMC algorithm suggest that the exponential autocorrelation 
times scale linearly (like a~^) if the length of the molecular-dynamics trajectories is 
scaled accordingly [19]. The algorithm however turns out to be non-renormalizable 
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in perturbation theory [6] and its scaling behaviour in the presence of interactions 
may therefore be completely different. 

The empirical studies reported later actually show that the HMC algorithm (on 
lattices with open boundary conditions) appears to fall into the universality class of 
the Langevin equation. In particular, the autocorrelation times scale approximately 
like rather than linearly. Prom this point of view, the non-renormalizability of 
the HMC algorithm in perturbation theory merely reflects the fact that the leading- 
order theory is in the wrong dynamical universality class and therefore not a suitable 
starting point for the perturbation expansion. 

3.5 Making QCD simulations safer 

In practice, numerical simulations should be much longer (by, say, a factor 100 at 
least) than the longest exponential autocorrelation time tq, as otherwise a proper 
sampling of the functional integral is not guaranteed and the simulation may conse- 
quently be biased in an unpredictable way. Usually the integrated autocorrelation 
times of the quantities of interest are monitored, but it should be noted that the cor- 
rectness of the simulation results ( within the estimated statistical errors ) cannot be 
taken for granted if only these autocorrelation times are much smaller than the total 
simulation time. 

Integrated autocorrelation times of both physical and other observables can in 
fact be very much smaller than tq. In particular, the autocorrelation times of noisy 
quantities (large Wilson loops, for example) tend to be practically unrelated to the 
exponential autocorrelation times. To illustrate this point, consider an observable 



where c is a constant and r/ a statistically independent Gaussian noise with mean zero 
and unit variance. Oq has the same expectation value as Oi and its autocorrelation 
function is given by 



At large c, i.e. when the added noise term is large, the integrated autocorrelation time 
'nnt(C'o) decreases like 1/c^ and can therefore be made arbitrarily small. Nothing 
is gained in this way, but the example shows that integrated autocorrelation times 
may not be representative of the true autocorrelations in the simulation. 

Exponential autocorrelation times are difficult to determine reliably if very long 
simulations are impractical. In this case, a pragmatic way to proceed is to look for 



Co = Oi + crj, 



(3.7) 



Tooit) = c'^Sto + Tnit). 



(3.8) 
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Fig. 1. Autocorrelation time of the density E at physical time L/2 in the SU(3) 
gauge theory, plotted as a function of the flow time t (cf. subsect. 4.3). The simulation 
data (points) were obtained on a lattice of size 32^ with spacing a = 0.05 fm and open 
boundary conditions, using the SMD0.3 algorithm. The line is a fit to the data of the 
form Tint = Co — cie"''^* with co — 94, while the leading exponential autocorrelation 
time in the even-parity sector is found to be about 100 in these simulations. 

observables with large integrated autocorrelation times and to take the maximum 
of the latter as an estimate of tq. The observables that provide the best probes for 
autocorrelations should be sensitive to the smooth modes of the gauge field, since 
these tend to be updated least efficiently. Moreover, for the reasons given above, 
good probes are likely to have small statistical fluctuations. Quantities obtained by 
integrating the Wilson flow [1], such as the average action density E at positive flow 
time, satisfy both criteria and are therefore recommended probes (see fig. 1). 



4. Numerical studies 



In order to verify and complement the theoretical discussion in the previous sections, 
we performed extensive simulations of the SU(3) gauge theory with open boundary 
conditions. The algorithms, observables and simulation parameters used in these 
studies are described in this section. 

4-1 Simulation algorithms 

Both the HMC and the SMD algorithm operate in phase space, i.e. on the gauge field 
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U{x,iJ,) and its canonical su(3)-valued momentum field Tr{x, fj,). The 0(a) boundary 
counterterm (2.14) is not included in the Hamilton function 



H{Tr,U)='^{7:,Tr) + SG{U) 



(4.1) 



of the system, partly for simplicity and partly because the term is unimportant in 

the present context. 

The HMC algorithm proceeds in cycles, where in each cycle one first chooses the 
momentum field randomly, with normal distribution, and then evolves the fields ac- 
cording to the molecular-dynamics equations that derive from the Hamilton function 
(4.1). In our simulations, the equations were integrated from molecular-dynamics 
time to r using tt-q iterations of the 4th-ordcr Omclyan-Mryglod-Folk (OMF) inte- 
grator defined through eqs. (63) and (71) in ref. [20]. At the end of the evolution, the 
fields are submitted to an acceptance-rejection step that corrects for the integration 
errors. This algorithm has two parameters, r and no, and requires the derivative of 
the gauge action to be calculated 5no times per cycle. 

In the case of the SMD algorithm, one proceeds in essentially the same way, but 
the momentum field is only partially refreshed according to 

Tr{x, /J.) -)• ci7r(x, n) + C2v{x, /i), (4.2) 



where f (x, //) is a randomly chosen momentum field with normal distribution, while 
7 and St are parameters of the algorithm. The molecular-dynamics equations are 
then integrated from to St by applying a single iteration of the 4th-order OMF 
integrator and the fields arc finally submitted to an acceptance-rejection step. When 
rejected, the fields are reset to their values before the integration, except for a change 
of sign 

7r(x,/x) ^ -7r(x,/x) (4.4) 

of the momentum field (see ref. [21] for a straightforward proof of the correctness of 
the algorithm). Note that the simulation time t elapsed after n SMD cycles is, by 
definition, equal to uSt, irrespectively of the rejection rate. 

Since the OMF integrator is applied only once, the molecular-dynamics evolution 
time St is usually set to a value much smaller than 1 in order to guarantee a high 
acceptance rate -Pace- Otherwise the SMD algorithm is frequently backtracking, on 



-jdr 



C2 = (l-c?)v^ 



(4.3) 
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average after every period of time equal to 

t.cc = 5T-^^, (4.5) 

and thus tends to become inefficient. With respect to the leapfrog and the 2nd-order 
OMF integrator, the 4th-order OMF integrator has the advantage that very high 
acceptance rates can be achieved with a moderate computational effort. 

4-2 Stochastic equation, parameter scaling and the SMD^ algorithm 

In the limit 6t — 0, the SMD algorithm amounts to solving the stochastic molecular- 
dynamics equations 

dsUs{x,ix) ='Ks{x,ix)Us{x,ii), (4.6) 

^s^^,{x,^x) = -T'^idl^ScKU,) - 2fioT^,{x, fi) + r],{x, /i), (4.7) 

where rjs a Gaussian random noise with mean zero and variance 

{v:{x,M{y, ^)) = 4.fioS'''6^J{s - r)a-Xy. (4.8) 

In these equations, the evolution time s and the mass /xq are related to the simulation 
time t and the parameter 7 through 

s = ta, Ho = 7/2a, (4-9) 

respectively. Evidently, eqs. (4.6), (4.7) reduce to the standard molecular-dynamics 
equations if fiQ is set to zero (see appendix A for the definition of the derivative of 
the gauge action). 

When the continuum limit is approached, the scaling behaviour of the simulation 
algorithms depends on how their parameters are scaled. The fact that the evolution 
time in eqs. (4.6), (4.7) has dimension [length] suggests to scale the HMC trajectory 
length r proportionally to l/a [19]. For the same reason, one can argue that fiQ 
should be scaled like a physical mass up to a logarithmically varying renormalization 
factor perhaps. This choice of the parameter scaling (which, however, leads to non- 
removable ultra-violet singularities in perturbation theory [6]) will be referred to as 
free-field scaling. 

Alternatively, if 7 is held fixed, and if 6t is such that the continuous-evolution time 
tacc is on the order of the exponential autocorrelation times (or larger), the SMD 
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algorithm effectively performs a numerical integration of the Langevin equation [6] . 
For clarity, we use the acronym SMD^ for the SMD algorithm with this parameter 
scaling. 

4-3 Observables 

As already noted in subsect. 3.5, observables based on the Wilson flow probe the slow 
modes of the gauge field and arc therefore well suited for studying autocorrelations 
in QCD simulations. A review of the Wilson flow is beyond the scope of this paper 
and we merely write down the differential equation 

dtVtix,fi) = -ag^oT'^id^^ScWtrnx,,,), T4(^,/x)|,=o = U{x,^^), (4.10) 

that generates the flow Vt{x,iJ,), t > 0, in the space of gauge fields (see refs. [1,4,22] 
for a comprehensive discussion of the flow and some of its surprising properties). 

In the course of the simulations, the observables are evaluated at fixed separations 
in simulation time. Starting from the current gauge-field configuration U{x,iJ,), we 
first integrate the flow equation (4.10) numerically up to some flow time t. The field 
tensor G^^(a;) of the gauge field Vt(x,//) generated in this way is defined through 
the clover formula, i.e. through the four plaquette Wilson loops in the (/x, ^')-plane 
that start and end at x (at the boundaries .xo = and a;o = T we set Gok(,x) = 0). 
The primary observables considered are then the time-slice averages 

^^^'^ = (4-11) 

X 

of the action density and the time-slice sums 

Q{xo) = Yl W-WGM^'(a^)Gpa(x)} (4.12) 

X 

of the topological charge density. Evidently, the autocorrelations of the total charge 

T 

Q = aJ2 (4.13) 

xo=0 

are studied as well. In all these equations, the dependence on the flow time has been 
suppressed for simplicity. 

At positive flow time t, the expectation values of arbitrary (finite) products of the 
observables E{xo), Q{xo) and Q do not require renormalization and are expected to 
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Table 1. Lattice parameters 



L/a 




a [fm] 


to/ a?* 


16 


5.96 


0.1000(6) 


2.698(3) 


20 


6.09 


0.0802(5) 


4.203(5) 


24 


6.21 


0.0667(5) 


6.086(7) 


32 


6.42 


0.0500(4) 


11.045(15) 


40 


6.59 


0.0402(3) 


17.49(4) 



* Calculated at physical time L/2 on the {L/a)'^ lattices 



have a well-defined limit when the lattice spacing is taken to zero [1,22]. While these 
expectation values do not have any obvious interpretation in terms of glueballs or 
colour flux tubes, for example, they are properties of the continuum theory which 
reflect the dynamics of the smooth modes of the gauge field (the smoothing radius 
being roughly equal to \/8i)- In particular, as explained in ref. [1], on lattices with 
periodic boundary conditions, the topological charge Q (as defined here) converges 
to an integer- valued observable in the continuum limit, which labels the topological 
sectors of field space. 

4-4 Lattice and simulation parameters 

In table 1 we list the spatial sizes and the inverse gauge couplings P = Q/g^ oi the 
lattices that we have simulated. The number of lattice points in the time direction 
(which is equal to T/a + 1) coincides with L/a in most cases, but lattices with larger 
time extent have been considered too. For the conversion to physical units, we use 
the Sommer radius ro = 0.5 fm [23] and the results obtained for ro/a by Necco and 
Sommer [24] . The values of the lattice spacing determined in this way (3rd column 
of table 1) are such that all lattices have the same physical size L, as is desirable for 
a scaling study. 

As a reference for the Wilson flow time t, we prefer to use the scale to determined 
through the implicit equation [1] 

{^'(^(^/2)>},=,„=0.3. (4.14) 

At flow time to, the Wilson flow has a smoothing range approximately equal to ro, 
i.e. this point in flow time is about where the non-pcrturbative regime sets in. Since 
L is quite small in physical units, the values of to/a'^ quoted in table 1 are probably 
affected by finite- volume effects and they are, in fact, a few percent lower than those 
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Table 2. Parameters of the HMC algorithm 



Lattice 


r 


no 


p 

^ acc 




JVcnfg 


16^ 


2.0 


6 


0.953 


6 


30697 


20^ 


2.5 


9 


0.975 


10 


25713 


24^ 


3.0 


12 


0.979 


15 


25625 




-i.O 


20 


0.985 


21 


21011 




Table 3. Parameters of the SMD0.3 


algorithm 




Lattice 


5t 




^acc 


At 


-^cnfg 


16^ 


0.1410 




516(2) 


5.92 


35093 


20^ 


0.1128 




748(2) 


9.14 


20209 


24^ 


0.0940 




1009(3) 


14.5 


20521 


48 X 24^ 


0.0818 




1205(2) 


13.7 


70000 


80 X 24^ 


0.0809 




964(2) 


13.6 


40991 


32^ 


0.0705 




1633(3) 


23.7 


20473 


40^ 


0.0561 




2.168(5) 


.18.6 


15976 



previously obtained in ref. [1] at L 2.4 fm. In the present context, the effect can 
however be safely ignored since L is the same on the lattices simulated. 

The trajectory length r in the HMC simulations was scaled according to the free- 
field parameter scaling, and the number no of integration steps (each consisting of 
one iteration of the 4th-order OMF integrator) was then tuned to achieve fairly 
high acceptance rates Pace (sec table 2). In a second set of simulations, we used the 
SMD-y algorithm. Some experimenting suggests that the autocorrelation times of the 
observables considered have a fiat minimum near 7 = 0.3 and we therefore decided 
to stick to this value of 7. The other parameter of the algorithm, 5r, was adjusted 
to ensure acceptance over an average simulation time face significantly larger than 
the exponential autocorrelation times (see table 3). 

The observables were measured at the separations A.t in simulation time quoted 
in tables 2 and 3. On each lattice, a fairly large number AT^nfg of configurations were 
analyzed, the total length of the simulations thus being equal to iVcnfgAt. 
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Fig. 2. Normalized autocorrelation functions of the observables E{L/2), Q{L/2)^ 
and at flow time to , plotted as a function of the simulation time lag t given in units 
of (L/af. The SMD0.3 algorithm was used all cases shown here. For better legibility, 
the data points obtained on the coarsest lattices (16^ and 20^) are coloured in grey, 
while the black points are those from the other lattices (24'', 32* and 40*). 

5. Simulation results 

5.1 Scaling properties of the autocorrelation functions 

While the chosen observables do not require renormalization, the flow time at which 
they are evaluated should be scaled like a physical quantity of dimension [length]^ 
in the continuum limit. In the following, the flow time is set to the reference time 
to, the results at other values of the flow time being similar as long as the short-time 
regime (where lattice effects are large) is avoided. 

The SMD0.3 algorithm is renormalizable to all orders of perturbation theory since 
it effectively integrates the Langevin equation (cf. subsects. 3.3 and 4.2). Moreover, 
with open boundary conditions, the topology barriers that otherwise slow down the 
algorithm are absent. It is therefore not unreasonable to expect that the normalized 
autocorrelation functions of the selected observables converge to universal functions 
in the continuum limit, provided the simulation time is scaled according to eq. (3.6). 

The autocorrelation functions plotted in fig. 2 in fact behave as expected if one 
assumes that the renormalization constant Zt varies only slightly on the lattices con- 
sidered. Note that all points obtained on a given lattice are statistically correlated. 
In particular, the seemingly systematic deviation of the measured autocorrelation 
functions on the 40^ lattice from those on the 32^ and 24"^ lattices may very well be a 
statistical fluctuation. Large deviations are however seen in the case of the time-slice 
and the total topological charge on the coarser lattices, where topology-tunneling 
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500 1000 1500 500 1000 1500 500 1000 1500 



Fig. 3. Integrated autocorrelation times of the observables E{L/2), Q{L/2)'^ and 
at flow time to, as obtained on the {L/a)'^ lattices using the HMC algorithm (open 
circles, scale factor Z = 1.32) and the SMD0.3 algorithm (full circles, Z = 1). Many 
HMC points lie on top of the SMD0.3 points and thus mask the latter. The curves are 
straight-line fits of the SMD0.3 data. 

transitions are not totally suppressed and thus reduce the autocorrelations. Langevin 
scaling then sets in as expected once these lattice artifacts become unimportant. 

On physically large lattices, the four-point autocorrelation function of the topo- 
logical charge Q is dominated by its disconnected parts. The normalized two-point 
autocorrelation function of is then related to the one of Q through 

PQ<t) {pQit)}' ■ (5.1) 

Although the simulated lattices are not very large in physical units, we found that 
eq. (5.1) is accurately satisfied. In particular, the autocorrelation functions of Q and 
scale in practically the same way. 

5.2 Autocorrelation times 

Similarly to the energy spectrum in finite volume, the exponential autocorrelation 
times depend on the symmetry sector considered. In particular, eq. (5.1) suggests 
that the longest autocorrelation time in the odd-parity sector is larger, by a factor 2 
perhaps, than the one in the even-parity sector. On the basis of the data shown in 
fig. 2, we estimate that the latter is about 1.2 x (ro/a)^ (thus ranging from 30 to 187) 
in the case of the SMD0.3 algorithm and the (L/a)^ lattices we have simulatedf. As 

f In accordance with the conventions adopted in sect. 3, all autocorrelation times are quoted in 
units of simulation time (i.e. molecular-dynamics time in lattice units). 
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Fig. 4. Integrated autocorrelation times of E{xq) and Q{xof at flow time to, plotted 
as a function of the physical time xq in lattice units. The data were obtained on the 
24"^, 48 X 24^ and 80 x 24^ lattices using the SMD0.3 algorithm. For better legibility, 
the data points obtained on the two smaller lattices are coloured in grey. 



usual, such estimates should be taken with a grain of salt, because the slowest modes 
in the system may not couple sufficiently strongly to the measured observables for 
their effects to be seen in the available data. 

The integrated autocorrelation times plotted in fig. 3 and their errors were cal- 
culated following the lines of appendix B. As is evident from the figure, the auto- 
correlation times all scale linearly in \/a? and thus as expected for algorithms that 
integrate the Langevin equation. Prom the point of view of the continuum limit, the 
intercepts at L/a = of the straight lines in fig. 3 are 0{a?) lattice corrections to 
the Langevin scaling, while the ratios of their slopes are universal properties of the 
simulation dynamics. 

Figure 3 also shows that the HMC algorithm (with free-field parameter scaling) 
scales like the SMD0.3 algorithm. The matching of the autocorrelation times requires 
a renormalization of the simulation time by the factor Z ~ 1.32, but in terms of 
computer time, HMC simulations tend to be faster than SMD0.3 simulations, because 
a very accurate integration of the molecular-dynamics equations is not needed. 

5.3 Dependence on the time-like extent of the lattice 

In practice, the time extent of the simulated lattices will often have to be larger 
than the one of the (L/a)^ lattices in our scaling studies. Autocorrelations in gen- 
eral depend on the physical situation and thus also on the lattice geometry. For 
illustration, the autocorrelation times of E and calculated on three lattices with 
the same spacing and spatial size, but different time extent T, are plotted in fig. 4. 
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Close to boundaries of the lattice, the autocorrelation times shown in these plots 
are thus practically independent of T, while well inside the lattices they rapidly con- 
verge to a constant value when T is increased. The behaviour of the autocorrelation 
times of these observables discussed in subsection 5.2 is therefore expected to be 
representative of the situation on larger lattices as well. 

The total topological charge Q is a special case, because it can only change (at 
small lattice spacings) by flowing in and out of the lattice. In the course of a simula- 
tion, the measured values of Q fluctuate around the origin with a standard deviation 
that increases proportionally to VtL^ on large lattices. The charge however flows 
through the boundaries with a rate proportional to \/Z^ only. The simulation time 
required for a significant change in Q must therefore be expected to grow with T 
(proportionally to T if Q performs a random walk). On the 24"*, 48 x 24"^ and 80 x 24'^ 
lattices, we actually find that the autocorrelation times of (42.1(2.5), 113(6) and 
148(10), respectively) grow roughly linearly with T. 

We wish to conclude this discussion by emphasizing that the autocorrelation times 
on lattices of a given physical size arc expected to scale linearly in 1/a?. Indepen- 
dently of the chosen geometry, the computational effort for HMC simulations with 
the standard leapfrog integrator, for example, thus scales approximately like 1/aJ . 



6. Conclusions 

The theoretical and empirical results presented in this paper show that the topology 
barriers in the SU(3) gauge theory can be avoided by choosing open boundary con- 
ditions in the time direction. Moreover, on lattices with these boundary conditions, 
the HMC and the SMD-y simulation algorithm both appear to fall in the dynamical 
universality class of the Langevin equation, i.e. simulations based on these algorithms 
slow down proportionally the square of the lattice spacing when the continuum limit 
is approached. 

In our numerical studies, the autocorrelation times of the topological charge (as 
well as those of observables unrelated to the latter) went up to values greater than 100 
in units of molecular-dynamics time. While such autocorrelations may be affordable 
in a given case, the experience suggests that there is ample room for algorithmic 
improvements. A separate treatment of the high-frequency and the smooth modes 
of the gauge field, for example, might be worth considering at this point. 

Open boundary conditions can easily be imposed in QCD with a non-zero number 
of sea quarks. We do not foresee any technical issues when these boundary conditions 
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are chosen, but an interesting theoretical question is whether the Langevin equation 
remains renormahzable in the presence of the pseudo-fermion fields that need to be 
introduced to be able to simulate the theory [27,28]. 

All simulations reported in this paper were performed on a dedicated PC cluster 
at CERN. We are grateful to the CERN management for funding this machine and 
to the CERN IT Department for technical support. 



Appendix A. Notational conventions 



The Lie algebra su(3) of SU(3) may be identified with the linear space of all traceless 
anti-hermitian 3x3 matrices. We choose the generators T", a = 1, . . . , 8, of the Lie 
algebra to be such that 

tr{r«r''} = -^S^K (A.l) 

The general clement X of su(3) is then given hy X = X"'T"' with real components X"" 
(repeated indices are automatically summed over). The Euclidean Dirac matrices 
7^, /X = 0, . . . , 3, are assumed to be hermitian. 

Gauge potentials Afj_{x) take values in su(3) and are normalized such that the field 
tensor and the covariant derivatives that appear in the Dirac operator are given by 

F^,. = d^A^ - d,A^ + [A^, A^l (A.2) 

D^ = d^ + A^. (A.3) 

On the lattice, the gauge-covariant forward and backward difference operators in 
presence of a lattice gauge field U{x,ii) act on the quark field ip{x) according to 

Vt,ip{x) = \{U [x, ii)ipix + afi) - ip{x)} , (A.4) 

V^ip{x) = i {ip{x) -U{x - ap,,ii)~'^ip{x - afi)} , (A.5) 
where a denotes the lattice spacing and fi the unit vector in direction /x. 
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The scalar product of any two vector fields oj{x, fj,) and v{x, ji) with values in su(3) 
is normalized such that 



{ijj,v) = — 2a^ tr{a;(x, /n)}. (A. 6) 



If J^{U) is a differentiable function of the gauge field, its derivative with respect to 
the link variable U{x,ijl) in the direction of the generator T" is defined by 



( e*'^''U{x,fi) if = (x,/^), 

t=o [ U {y, u) otherwise. 

(A.7) 



In particular, in the case of a scalar function T{U), the combination T"'d^^^T{U) is 
a vector field with values in su(3) that transforms under the adjoint representation 
of the gauge group. 



Appendix B. Calculation of integrated autocorrelation times 

The integrated autocorrelation times of the selected observables Oi were obtained 
as usual from the empirical estimates 

_ Ay- *tot-t _ _ 

TiiW = T — 7 yz (^^(^) - ^0 i^ii^ + - (B.i) 

of the autocorrelation functions Tii{t), where ttot = -^cnfgAt denotes the total simu- 
lation time of the run and Oi the average of the measured values of Oi. In all cases, 
the autocorrelation functions are found to decay exponentially at large time sepa- 
rations with remarkably consistent values of the exponential autocorrelation times. 
The estimate 

rint(Oi) ^ iAi + AtY, PiikAt), pi{t) = (B.2) 
k=i ^"y^> 

therefore rapidly approaches a constant value when the "summation window" W = 
kraa.xAt IS Sufficiently large. 



23 



On the {L/aY lattices considered, the summation window for even-parity observ- 
ables was set to 

6.0 (HMC runs), 

(B.3) 

4.5 (SMDo.3 runs). 

Given the measured exponential autocorrelation times (subscct. 5.2), the systematic 
error that derives from the truncation of the sum (B.2) is estimated to be at most 
3% with this choice. The statistical errors of the autocorrelation functions and the 
integrated autocorrelation times were determined using the Madras-Sokal approxi- 
mation [25] (see ref. [26], appendix E, for a detailed description of the procedure). 



W = ro/a 
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